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Abstract 

We present a new formulation of the correlated electron-ion dynamics (CEID) scheme, which sys- 
tematically improves Ehrenfest dynamics by including quantum fluctuations around the mean-field 
atomic trajectories. We show that the method can simulate models of non-adiabatic electronic tran- 
sitions, and test it against exact integration of the time-dependent Schrodinger equation. Unlike 
previous formulations of CEID, the accuracy of this scheme depends on a single tunable parameter 
which sets the level of atomic fluctuations included. The convergence to the exact dynamics by 
increasing the tunable parameter is demonstrated for a model two level system. This algorithm 
provides a smooth description of the non-adiabatic electronic transitions which satisfies the kine- 
matic constraints (energy and momentum conservation) and preserves quantum coherence. The 
applicability of this algorithm to more complex atomic systems is discussed. 

PACS numbers: 31.15.Qg, 31.50.Gh 
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I. INTRODUCTION 



Ordinary molecular dynamics (MD)^ uses Hamiltonian equations of motion (EOM) to 
describe the evolution of an atomic system. The validity of that approach relies on the 
Born-Oppenheimer (BO) approximation which states that — in most cases — the atomic 
motion induces a slow (adiabatic) perturbation of the electronic dynamics. Therefore, if 
the system is originally prepared in an electronic eigenstate (e.g. the ground-state) for a 
given atomic configuration, it will evolve into the same electronic eigenstate for the evolved 
atomic configuration. Furthermore, according to the BO approximation, the ground-state 
electronic energy for a given atomic configuration can be employed as an effective (i.e. low 
energy) potential energy surface (PES) for the atomic motion. 

The BO approximation is not longer applicable whenever electronic transitions between 
surfaces are relevant: for instance, where there is a non-radiative transition following an 
earlier photo-excitation. However, it is still possible to extend the scope of MD by allowing 
more than a single PES — one for every electronic state — and by providing a meaningful 
way to make transitions among PES i.e. non-adiabatic electronic transitions. 

The extension of ordinary MD to deal with processes involving many PES has been exten- 
sively pursued over recent decades and some effective algorithms are available for atomistic 
simulations. The most used are: Ehrenfest dynamics (ED),^'^ molecular dynamics with 
quantum transitions (MDQT),^-'', mixed quantum-classical dynamics (MQCD),-'^ and ab 
initio multiple spawning (AIMS)^'^. They have been employed to simulate quantum dissi- 
pative dynamics^'^° as well as non-adiabatic processes such as photo-chemical reactions,— 
polaron formation in conjugated polymers,—, and proton transfer in solutions^. Although 
less successful in practical applications, we also mention other dynamical schemes that pro- 
vide valuable theoretical insights: the frozen Gaussian approximation,^^ which is one of the 
pillars of AIMS, and the Gauss-Hermite wave-packet expansion,— which shares similarities 
to the method presented in this paper. 

In order to extend the scope and the accuracy of the aforementioned methods, a new 
approach called correlated electron-ion dynamics (CEID)^^^, has been introduced. This 
method has been mainly applied to study the heat production and dissipation in model 
metallic nanostructures,-ii^»i^ a problem relevant for nanotechnology. Indeed, other algo- 
rithms based on non-smooth dynamics, that is, those which allow for either sudden surface 
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hopping (like MDQT) or sudden wave-packet spawning (like AIMS), are expected to be 
less efficient in the simulation of systems with a dense, gapless electronic spectrum, like 
a metal. That is a consequence of the slight adjustments of the average atomic positions 
and/or velocities that might need to be imposed after a sudden transition in order to ensure 
total energy and momentum are conserved.-'^ If the number of crossings is large, the com- 
putational overhead due to these adjustments might be non- negligible. MQCD, although 
it is often implemented by using both surface hopping and mean-field evolutions, is based 
on a sound theory^i^i and conserves the kinematic constraints. On the other hand, time- 
translation invariance is valid only approximately in MQCD^ and numerical instabilities 
have so far limited the surface-hopping implementation of MQCD to relatively short time 
simulations.™ Finally, ED — although it evolves smoothly — has been shown to poorly 
describe the atomic heating caused by electron-ion interaction.— 

A further complication arises when a quantum sub-system is coupled with large quantum 
reservoirs (which are not treated in detail), as for a nanostructure connected to macroscopic 
leads. In this case, it is hard to define meaningful PES in terms of the sub-system degrees 
of freedom only, and so algorithms based on the surface hopping paradigm are expected to 
be less accurate. This problem is absent if the quantum sub-system is coupled to a classical 
dissipative environment; in this case MQCD or its latest variant^^ can give reliable results. 

In principle, the CEID EOM form an exact, yet infinite, kinetic hierarchy which corrects 
ED by means of the so-called small amplitude moment expansion (SAME) of the Liouville 
equation.- So far, a few schemes to truncate the hierarchy have been proposed^ in order 
to simulate the CEID EOM and currently available algorithms are restricted to a mean- 
field second moment approximation.— Although those algorithms are accurate enough to 
describe — at least qualitatively — nanostructure heating, the existence of a practical trun- 
cation scheme which converges to the exact quantum dynamics has not been demonstrated 
until now: in this paper we show that a convergent truncation scheme actually exists and 
we provide a new practical CEID algorithm whose accuracy depends on a single tunable 
parameter. We provide a validation of our theory by comparing CEID results against exact 
integration of the time- dependent Schrodinger equation for a model two level system (2LS) 
— i.e. the simplest system which displays non-adiabatic transit ions^^i^^ — in two contrasting 
parameter regimes. 

The rest of the paper is organized as follows: In Sec. Ull we introduce the physics of 
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the model 2LS employed to test our new CEID algorithm. In Sec. Illll we describe the exact 
algorithm by which we produced the benchmark calculations, while in Sec. llVl the new CEID 
scheme is derived in detail. Finally, numerical results from simulations of the model 2LS 
are collected in Sec. |V] and the conclusions and perspectives of this work are discussed in 
Sec. ED 

II. A CASE STUDY: THE TWO LEVEL SYSTEM 

In this section we introduce the model 2LS that we employed to test the convergence 
properties of our CEID scheme. [Numerical findings are reported in Sec. |Vl] Here we 
discuss the physics we expect to address before explaining the details of the algorithms we 
use. 

It is widely recognized that a one- dimensional 2LS illustrates many of the fundamental 
features of a non-adiabatic system and, at the same time, it retains the simplicity of a low- 
dimensional model.— 1^ In general, the Hamiltonian for a system made by electrons and ions 
(in the absence of external fields) can be written as follows: 



where R and P are the quantum operators for the atomic position and momentum, while 
the electronic dependence of H is collected into Hf,. In particular, the first term in the RHS 
of Eq. ([1]) accounts for the atomic kinetic operator while a sort of atomic potential operator 
is described by the second term. 

There is a lot of freedom in constructing the potential term, ifg? but in this paper we 
focus only on the following parametrization: 



which describes two parabolic PES (diagonal entries) linearly coupled through a kind of dipo- 
lar interaction (off-diagonal entries). This is a non-adiabatic representation of the electronic 
PES in which the electronic basis is independent of the atomic coordinate. The adiabatic 
representation can be obtained as usual by diagonalizing He{R). Eq. ([T]) and ([2]) depend on 
the following parameters: the atomic mass M, the harmonic constant K, the electron-ion 
coupling constant fc, the surface displacement Rq, and the PES energetic offset As. 



H = P^/2M + He{R) , 



(1) 
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Since both the PES are confining, we expect to see periodic electronic transitions between 
the two PES driven by the electron-ion interaction. For instance, a state can be prepared 
as the atomic ground-state on the upper electronic PES. The atomic position R experiences 
quantum oscillations and so the system cannot be exactly localized in the minimum of the 
PES {R = according to our parametrization) . Since this initial state is not an eigenstate 
of the interacting Hamiltonian (i.e. for /c 7^ 0), the system will eventually make a transition 
into the lower PES. We stress that this process must conserve the total energy so that an 
atomic transition must accompany the electronic transition. For instance, the electronic 
process described above can be viewed as an initial decay since the atomic potential energy 
is effectively decreased. Therefore, an increase of the atomic kinetic energy is expected as a 
consequence of this decay in order to conserve the total energy. This, in a nutshell, is the 
heating of an atomic degree of freedom caused by the electron-ion interaction.— However, 
we stress that the aim of the present work is not the study of a quantum decay process, but 
the illustration of the convergence properties of a new CEID algorithm. Other models must 
be used to address the physics of quantum thermalization. For instance, the spin-boson 
model^ describes a 2LS coupled to an environment made by a collections of many quantum 
harmonic oscillators. This model can be effectively simulated^i^ and it might provide a 
future test case for our new CEID algorithm. 

There is an interesting general feature displayed by the electronic dynamics of our model 
2LS. According to the initial condition described above, the electronic transition is due to 
the quantum fluctuations only, because the dipolar interaction is exactly zero for a classical 
atom perfectly localized in the minimum of the upper PES. On the other hand, quantum 
fluctuations are completely neglected in the sort of mean-fleld description of the atomic 
motion employed in ED. As a consequence, we do not expect ED to reproduce the initial 
transition from the upper PES to the lower PES, which can be thought of as spontaneous 
phonon emission. If conflrmed, this behavior will provide further evidence that the exchange 
of energy from the electronic to the atomic degrees of freedom cannot be properly addressed 
by ED.-^ 

In order to be as transparent as possible when discussing the dynamical features of our 
model 2LS, we choose to measure the values of the parameters in the Hamiltonian in terms 
of natural units. The natural energy scale is given by the harmonic quantum, hu, where 
uj = ^yK/M, and so the time can be measured in units of the harmonic period, 271 /uj. 
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Introducing a mass scale, Mq, (its actual value is not important here), a length scale and 
a linear momentum scale are immediately obtained: ao = \J h/ {Mquj) and bo = y/MoHuo, 
respectively. 

For our purposes, not all the parameters need to be varied during the numerical exper- 
iments. First of all, we fixed both the atomic mass (M = Mq) and the harmonic constant 
{K = Mquj"^) and then we took the PES offset to be equal to one harmonic quantum 
{Ae = tujj). This is equivalent to saying that the atomic ground-state on the upper PES, 
IXo"'*)) has got exactly the same energy as the first harmonic excitation on the lower PES, 
IXi'"*)- [Here we are neglecting the electron-ion coupling and using the notation introduced in 
appendix [Dl] If the coupling constant is small enough (see appendix [Dl) , the time-evolution 
of our 2LS can be understood starting from these two states. We confined this study to 
this weak coupling regime, i.e. /c < 0.1 Tiw/aQ. Finally, we report in this paper the numer- 
ical results for two different 2LS geometries only: the unshifted 2LS, with Rq = 0, and the 
shifted 2LS, with Rq = cto- The PES for these two cases are shown in Fig. [l](a) and Fig.[2](a), 
respectively. Although many other 2LS geometries have been studied, e.g. with larger Rq 
or with other values of K (including the possibility of different harmonic force constants for 
the lower and the upper PES), they gave numerical outcomes qualitatively similar to either 
the shifted or unshifted 2LS, and so the details are not reported here. 

A linear combination of the two low-lying resonant states might be employed to give a 
qualitative account of our model 2LS dynamics i.e. the wave-function at time t might be 
approximated as: 

m ^ comi""^) + c,{t)\x?) , (3) 

where cq and ci are time- dependent complex coefficients. This state would give a distri- 
bution of the atomic position R more or less localized around each of the two classical 
equilibrium positions, namely R = for the upper PES and R = Rq for the lower PES. 
Eq. ([3]) might not be a good ansatz for the evolved state as soon as the displacement Rq is 
large enough. Indeed, even a classical atom can be found far from its equilibrium position 
whenever this is energetically allowed. In this case, we should be able to describe an atomic 
wave-packet almost localized far from both i? = and R = Rq. Therefore, we will produce 
a manifest physical inconsistency if we assume that the generic 2LS state can be always well 
approximated by Eq. ([3]). This physical inconsistency can be easily fixed by taking a longer 
expansion of the exact evolved state ijj{t) in terms of the harmonic excitations of the atomic 
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degrees of freedom. Unfortunately, there is a cost to pay: the longer — and so the more 
accurate — the expansion, the more time-consuming will be the simulation. [See Sec. IVII 
for a further discussion of this point.] 

In the next two sections, two different ways to compute non-adiabatic dynamics are 
considered in detail. In Sec. IIIII a method based on the numerical diagonalization of the 
Hamiltonian of the full system is presented, while a new formulation of CEID is introduced 
in Sec. HVl 



III. EXACT INTEGRATION 



The 'exact' method of integration is based on numerical diagonalization of the full Hamil- 
tonian matrix and subsequent time evolution exploiting the eigenvectors and eigenvalues ob- 
tained in the diagonalization step. While it is highly accurate, in practice, such an approach 
naturally has to be limited to systems with a small number of degrees of freedom. However, 
for our purposes it provides an ideal scheme for producing benchmark results. A central 
quantity in this paper is the Wigner transform (WT) of an operator. Originally, the WT of 
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as 



a wave-function ip was defined in Ref. 

W{R, P) = J r{R + s)iiR - s) exp ■ s)d"s 

where R = . . . , P = (Pi, . . . , P„), and the time-dependence has been suppressed. 
The latter form in Eq. (jlj) is more common nowadays. In a straightforward extension of this 
definition, for an operator A, expandable in basis states {(pa}, A = J2ab^a{x)A°'''(f)l{x), the 
WT is given by 

A^{R,P) = ^A'^'^^ j <Pl{R + \s)<Pa{R - ^^)expQp ■ sy-s. (5) 

Of particular interest to us here, however, are WT with respect to the degrees of freedom of 
a subsystem, i.e. partial WT. These appear naturally when the system can be divided into 
subsystems on physical grounds. In this paper, where the system is a molecule, the system 
can be split into an electronic subsystem and the subsystem of the ions/nuclei. Consider 
an operator C acting on a system divisible into subsystems with basis sets and 
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The operator is assumed to be expandable as C — X^ABab 1^) ® 1^) C-^^"'^ {h\ ® {B\, and its 
partial WT with respect to the |A)-subsystem is 

(7^,^(7?, P) = ^|a)Cj^(it:,P)(6|, (6) 

ab 

where, with the position representation of 1^4), 

Note the important fact that Cw,a is still an operator in the |a)-subsystem. As in the rest 
of the paper the only WT used are partial WT with respect to 'ionic' or 'atomic' degrees of 
freedom, we shall henceforth write (7^ instead of Cyj^A, and also refer to the |74)-subsystem 
as the atomic and the |a)-subsystem as the electronic subsystems. 

The time-evolution of the system is generated by a Hamiltonian, whose eigenvalues and 
eigenvectors shall be £n and |^n), respectively. The basis states {|a)} and of the 

subsystems are taken to be time-independent from now on. We can expand an eigenstate in 
the product basis \^n) = '^Aa^n"' k) ® 1^)) or vice versa, \a) ® \A) = J2n'^Aa l^n)) where 
wc have the relations X^Aa^n^-^Xi — Ylm^n'^'^^b — ^bI- The expansion coefficients 

C and T> are time-dependent, 

CW=C„^"(0)exp(-^£„t); (8) 

in any particular case the numerical diagonalization will provide us with the £n and the 
coefficients C^"(0), which make up the eigenvectors, i.e. C^"(0) is the Aa-component of 
eigenvector n. 

We now can expand an operator G in the eigenstates or the product states: 

G = Y,\'^n) G^^ {^>m\ = 5^ |a) 1^) G^'"'\t) {h\ ® {B\ , (9) 

mn ABab 

with 

^ABab^^^ = (0)G"-(C)*(0)exp {-^{Sn- Sm)t) . (10) 

nm 

The partial WT with respect to the atomic subsystem is then 

G^{R, P,t) = Y^ l«) (i?, P, t) {b\ , (11) 

ab 
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where in turn 

G^^{R,P,t) = Y,^BA{R,P)G^'"'\t). (12) 

AB 

and 

Specializing to the particular case of the 2LS and its actual implementation on a com- 
puter, we introduce dimensionless quantities using the scale factors mentioned in Sec. [Ill In 
particular we obtain the dimensionless wave-function v^(0 = \/o^^{C'oC)i the dimensionless 
version of JF, Fnm{^,v) = ^^nm{o,o^,bof]) S'lid dimensionless eigenvalues En = Enlihjj). Us- 
ing Eq. flTU]) . the factoring Eq. f[T^ . and Eq. f[T^ . all the components of can be computed 
as functions of R, P, t (or their dimensionless counterparts). Note that this does not involve 
the numerical solution of a (potentially partial) differential equation, so we are not required 
to advance over many small time-steps in order to reach a given value of t. Limitations are, 
however, introduced by the need to truncate the oscillator basis to a finite number of states. 
As the purpose of the exact approach is to provide a benchmark for the CEID method, 
some care has to be taken to avoid truncation errors here. First one has to decide how 
many product and eigenbasis states to use in expansions like Eq. flTUl) . say A^. Then, the 
diagonalization has to be carried out using a number of product basis states M > N such 
that the lowest eigenvalues and corresponding eigenvectors are well converged; typically, 
we used A^ ~ 120 and M ^ 5N. The operator one is considering (in what follows it will be 
the density operator) should only have negligible coupling between states with index equal 
to or less than A^' (A^' < A^) to states with index larger than A^'. These first A^' states, 
which in the diagonalization step are produced as M-component quantities, should not have 
significant contributions from product basis components with index greater than A^. Only 
if these conditions are met, and the initial conditions for the operator are chosen to involve 
only the first A^' states, can we consider the numerical results for the time-evolution reliable 
and an adequate benchmark for CEID. The quality of the choice of A^ can therefore only be 
assessed after diagonalization. 

From the results produced for the purpose of comparison we show the occupations Na, 
a G {1,2}, of the electronic levels, and expectation values of position, momentum and the 
variance of the position, as functions of time. These have been calculated via the WT 
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p^{R, P, t) of the density operator p of the system, by numerical evaluation of the integral 

Nait) = jj pl^R, P, t)dRdP, (14) 
in case of the occupations, and of 

2 

(/(i?,P,t)) = J2 ff fiR,P,t)pl^{R,P,t)dRdP (15) 

a=l ^ ^ 

for f{R,P,t) = R, P, (i? — {R){t)Y, respectively, in the rest of the cases. 



IV. CORRELATED ELECTRON-ION DYNAMICS 

In this section we describe a new formulation of the correlated electron-ion dynamics 
(CEID) while the original formalism can be found in Ref. fiol. [See also Ref. 17 and Ref. [s] 
for further details.] For the sake of simplicity, the new CEID EOM has been derived here 
only for the one- dimensional case although multi-dimensional EOM are also known.— 

We start from the well known quantum Liouville equation: 



H,p 



(16) 



which is the EOM for the density matrix p of the system. [We use a dot to indicate time- 
derivative.] Unfortunately, a direct integration of Eq. (fT6|l is exceedingly time-consuming 
because it scales approximately as the cube of the Hilbert space dimension which is very 
large in most cases of interest. On the other hand, since atoms are much heavier than 
electrons, an expansion of their motion around the classical trajectories is often justified. 
It turns out that this kind of expansion cuts off the quantum fluctuations of the atomic 
degrees of freedom and so it effectively reduces the Hilbert space dimension. For instance, 
simulations of the semi-classical limit of Eq. f|T6l) ^^ have been shown to reproduce — at least 
qualitatively — the correct non-adiabatic dynamics of a few interesting test-cases. 

More generally, the density matrix p can be partially expanded with respect to the atomic 
degrees of freedom by means of a complete orthonormal system (COS).- By using the stan- 
dard Dirac's bra and ket notation, this expansion can be expressed as: 



oo oo 



P = X] 5Z I'l^n) Pn,m{<Pm\ , (17) 



n=0 m=0 
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where the functions {(pn} are a COS in the atomic subspace. As a consequence, in Eq. (fT7|) 
the electronic degrees of freedom are included in the matrix coefficients pn,m- A natural 
choice dictated by the kind of 2LS physics introduced in Sec. [TTl (i-e. confining PES) is to 
use the simple harmonic oscillator (SHO) eigenfunctions as atomic COS. We stress here 
that, although these functions are usually centered around a classical equilibrium point, any 
other reference point can be taken instead (see below). 

Since all the observables can be expanded as in Eq. f|T7|) . all the operations involving 
observables (e.g. averages) can be worked out by means of the observable matrix coefficients 
only. [In general, the matrix coefficients of the observable A are given by: An^m = {4>n\A\^rn) ■] 
For instance, the total energy of the system is given by: 

oo oo 

Etot = Tr {^p} = J] J] Tre [Hm,nPn,m} ■ (18) 
n=0 m=0 

[The two traces, Tr and Tre, apply to different linear spaces, namely the whole Hilbert space 
and the electronic subspace: see Sec. IIIII ] As a further example, the EOM for pn,m are 
obtained by plugging Eq. ( ITTl) into Eq. ( fT6l) : 

^ oo 

Pn,m ^ ^ |^-^n,A; Pk,m Pn,k Hf^ ^yi . (-^^) 

k=0 

It is easy to prove that Etot is in fact a constant of motion by taking the time-derivative of 
Eq. f|T8|) and then by using Eq. f|T9|) . 

The complete set of EOM, Eq. (fT9l) . cannot be directly simulated because it is not finite. 
Therefore, we must make an approximation and, quite naturally, we set to zero (as a matrix) 
every matrix coefficient with indices greater than A^, the CEID order. Nevertheless, after 
this truncation, the EOM are still fully quantized (and expressed by a proper Lie bracket), 
but they are restricted to a smaller Hilbert subspace. 

As for the exact scheme described in Sec. Illlt in order to make the classical limit of 
Eq. ( fT9l) more manifest, we use a partial Wigner transform (WT) i.e. a WT taken only 
with respect to the atomic degrees of freedom. Therefore, the partial WT of the operator A, 
Aw{R, P), is still an operator in the electronic subspace but it explicitly depends on what are 
now the classical atomic position R and momentum P. In the context of non-adiabatic MD, 
similar partial WT have been already consideredf^"^ and they have been shown to provide 
correct numerical results. In order to avoid confusion, we stress that, although WT seems 
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to map a quantum operator into a classical distribution, the dynamics remains non-classical 
because the WT of the product of two operators is in general not the product of the WT of 
the two operators: 

{Ab)^ = A^^B^ , (20) 
where is the non-commutative Moyal product^'^ 



exp 



^-^ [ d B,d p ~ d pd B, 



(21) 



[The arrows indicate the directions in which the derivative operators act.] The WT of any 
operator can be expanded as in Eq. f|T7|l by using the WT of the basis operators |0„)(0m|- 
As a consequence, the matrix coefficients of an operator are the same both in the original 
Hilbert space and the in transformed one. 

The WT of the Liouville equation can be formally stated as: 

Pw{t) = (Hw * Ptu - P«; * hJ^ . (22) 

It can be shown that from Eq. ( !22ll the usual Hamilton-Ehrenfest equations (i.e. the EOM 
for R = Tr{i?p} and P = Tr{Pp}) can be derived:- 

R = P/M , 

r / - \ ^ (23) 
P = F = -Tr{(^)p}. 

It is also reasonable to take the phase-space trajectory {R{t), P{t)) as a zero-order approx- 
imation of the true atomic dynamics if the quantum fluctuations are not too large. [Th e 



semi-classical limit of the WT is explained more extensively in Ref. [sil (see also Ref. 34).] 
This fact suggests that instead of taking the origin as a reference point in the phase-space 
{R, P) (i.e a fixed reference frame), one can take advantage of Eq. (1231) and use {R(t), P(t)) 
as reference (i.e. a mobile reference frame). After this mobile reference frame transform, 
the EOM becomes: 



in 



Hyj 'k Pm p^ -k Hy^ 



dpw \ P , f dp 



F . (24) 



OR J M \dP 

Although it is not apparent from Eq. fl2^ . the EOM can be still expressed by a proper Lie 
bracket — as in Eq. (122|) — by means of a different time-translation generator i.e. by a 
different Hamiltonian. [Mathematical details can be found in appendix [Bl] 
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At this stage two paths can be followed. In the first case, the Moyal product is expanded 
in h (usually up to first order^'^^) and a quantum-classical extension of the Liouville equation 
is obtained. However, although this approach is physically appealing, one must be aware that 
the EOM obtained this way cannot be formulated through a proper Lie bracket,—!^ and a 
generalized non-Hamiltonian bracket should be introduced instead. As a consequence, the 
evolution of a composite operator, [y4S](i(:), might be different from the composition of the 
separated evolutions, A[t)B[t),—^ (because the non-Hamiltonian bracket does not define a 
proper derivative), and the conservation of dynamical symmetries might be a problem^ (due 
to the violation of the Jacobi identity). It must be also recalled that the difference between 
the non-Hamiltonian and Lie brackets is only of order (9(/i)^° and that the non-Hamiltonian 
structure arises in a quite natural way for open systems, whether classical or quantum.— 

By following the other route, one uses the exact expression for the Moyal product and 
takes the WT of Eq. (ITTI) as a natural way to truncate pw [We have found a direct expansion 
in the transformed space — e.g. by using weighted orthogonal polynomials — to cause 
dangerous instabilities in the truncated dynamics.] Therefore, the action of the truncation 
super-operator Tyj can be defined as follows: 

[p^R, P, t)] = [p^{R + AR,P + AP, t)] 

N N 

= Y.Y1 ^' t)Pn,m{^R, AP) , (25) 

n=0 m=0 

where Pn,m{AR,AP) is the WT of \4>n){4'm\ in the mobile reference frame. [Properties of 
these functions are given in appendix |Al] We opted for this approach because it is still a 
full quantum scheme — but in a truncated Hilbert space — and the EOM for the density 
matrix can still be formulated by means of a proper Lie bracket (see appendix [Bl). 

Although an analytical expression of Pn,m is known, it is far more convenient to state 
a set of recurrence relations by taking advantage of the well-known SHO algebra. Details 
can be found in appendix (Bj here we report only the final form of the CEID EOM for the 
matrix coefficients pn,m (in the mobile reference frame): 
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Pn,m = [V + 2)(n + l)pn+2,m " (2/2 + l)pn,m + \/n{n - l)p„_2,m + 

-^Jm{m - l)pn,m-2 + (2m + l)p„,m - A/(m + 2)(m + l)p„,m+2j + 

+ ^ (i?) ,p„,„] - ^ ^AF (i?) y^±ip„+, „ + AF [R] y|p„_i^^+ 

-^Pn,m-iAF (R) - ^^^^Pn,m+iAF (i?) ^ + ^ (^^T (i?) ^(n + 2) (n + l)p„+2,m+ 

+i: (i?) (2n + l)p„,^ + k [R] — l)Pn-2,m — a/ m(m — l)p 

-(2m + l)p„,„i: (i?) - ^/{m + 2){m + l)pn,m+2k (i?)) , 

(26) 

where F = —dHe/dR, AF = F — F, and K = d'^He/dR^. [Terms involving higher deriva- 
tives of He should also appear in Eq. ( |26l) . but vanish in this case since the 2LS Hamiltonian 
we want to study is quadratic — see Eq. ([2]).] We recall that, according to our truncation 
scheme, one must neglect in the RHS of Eq. fl26|) those matrix coefficients whose indices are 
greater than the CEID order. Those equations, along with Eq. (!23|) . have been used to sim- 
ulate the 2LS dynamics described in Sec. [TTl In particular, the current implementation uses 
a second order Runge-Kutta non-adaptive algorithm to integrate Eq. and the standard 
velocity-Verlet algorithm to integrate Eq. ( 123|) . [A time-step At = 10~^/27r in our natural 
unit (see Sec. [TTl) been found to be appropriate for the precision required by the compar- 
ison between CEID an exact approaches reported in Sec. |Vl] At every integration step, the 
averaged coordinates, R and P, are evolved according to Eq. fl23l) for half a time-step, then 
the matrix coefficients are propagated through Eq. (12^ for a whole time-step, and finally the 
averaged coordinates are evolved by another half time-step. We verified that the accuracy 
achieved by this kind of symmetric Trotter decomposition is greater than what is obtained 
by means of a single evolution of the averaged coordinates for a whole time-step followed 
(or preceded) by a matrix coefficients propagation for a whole time-step. Numerical results 
can be found in Sec. |Vl 

We now briefly discuss the link between our new formulation and the original CEID. 
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A. Comparison with former CEID integration schemes 



At variance with the scheme described so far, the original formulation of CEID makes 
use of a completely different expansion which directly provides EOM for the moments of the 
density matrix.-ii^ The most relevant CEID moments are: Pe = Tr^ {p}, pi = 
and Ai = Tr^ |APp|, where Tr^ is the partial trace with respect to the atomic degrees of 
freedom, Ai? = R — R, and AP = P — P. [Higher order moments must be carefully defined 
because Af^ and AP do not commute.] On the other hand, analogous objects can be also 
introduced in the new formulation: 

f^n,m{t) = ^ / d^dP AP"AP"^p^(P, P, t) . (27) 

By using the property of the WT, it is easy to find a link between the new and the original 
notation: /io,o = Pe, Pq,i = Ai, and pifi = pi. Similar relations for higher CEID moments 
can be stated, but some extra attention must be payed in the derivation due to the non- 
trivial commutation relations between positions and momenta. It is worth noting that 
CEID moments provide valuable information about the system. For instance, the quantities 
(po,o)n,n = (Pe)n.n givc the probability of observing the system on the n-th PES and the 
average force (see Eq. fl23|) ) can be easily computed from F = Tr{F{R)pQfi} — TT{K{R)pi o}. 
Higher moments can be used to study electron-ion correlations.— 

Moments defined in Eq. ( l27l) can be expressed in terms of the matrix coefficients by means 
of the following linear transform: 

Pn,m = J2A^;rPr,s, (28) 
r,s 

where 

Apf = [ dPdPAP"AP"^P,,(P,P) . (29) 
zTin J 

As usual, a set of recurrence relations for A"^'^ can be found and — at least in theory — 
CEID moments of any order can be computed. In practice, only the low lying moments are 
relevant and here we give a short selection of them: 
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P2fl = +y X [Vn(n - l)p„,„_2 + (2n + l)pn,n + \/ n{n - l)pn-2,n 



(30a) 
(30b) 
(30c) 
(30d) 
(30e) 
(30f) 



n=0 



As anticipated in Sec. [H the zero order CEID (i.e. for = 0) is equivalent to the ED:— 

1 



P0,0 — -r 

in 



He (R) , po,o 



(31) 



[We have used Eq. ( l30l) (a) and Eq. ( l26i) .] We also stress that, in this case, /io,i = Pi,o = 
Pi^i = and that both /io,2 and /i2,o are proportional to /io,o- 

To clarify the link between the new and the original formalism, it is helpful to write down 
the first order (A^ = 1) EOM in terms of the CEID moments. This can be done by inverting 
Eqs. ( I30|) (a-c.f) to express po,o, Po,i, Pifl, and pi^i as functions of po,o; Po,i) Pi,0; and /t2,o- 
[At the same CEID order pi^i = and that /io,2 = (^o/'^o)A2,o-] The final result is: 
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-\{AF{R),fi2,o 
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(32b) 



Pi,o 



1 1 

— P0,1 + -T 
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F{R),fio,^ 
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(32c) 
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ih 



He{R),fl2,l 



+ 



+ 



+ 



2bl 



AF(/2),/io,i 



(32d) 



These equations might be compared with the Eq. (8) of Ref. Il7| (the mean-field second 
moment approximation) keeping in mind that in that paper the following ansatz as been 
made: /i2,o = C^'^fio^, fii^i = C^'^fio^o, and /io,2 = C^'^P'Ofi, where C^'^, C^'^, and C^'^ 
are time-dependent quantities. According to our initial conditions (see Sec. [Tll), the initial 
values of these variables are: Cr^r^O) = al/2, Cr^p{0) = and Cp^p{0) = bl/2. In order 
to find the EOM for Cr^r, one can trace Eq. (1321) (d) and it turns out that, to the first 
order in the coupling constant fc, Cr,r = 0. Remarkably, by assuming that the matrix K is 
proportional to the unit matrix and by substituting /i2,o = (cto/2)/^o,o into Eq. fl32p (a-c). we 
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obtain EOM for /io,o, P'0,i, and /ii_o which are equal to the ones stated in Eq. (8) of Ref. 
up to the first order in the coupling constant. Although there is no reason to believe that 
this agreement must be restricted only to a given set of initial conditions, it is not clear yet 
how it might be proved right for higher CEID order or different basis set expansion. 



B. Energy conservation 

By using the original CEID scheme, it is possible to write the total energy, Eq. (ITS!) . in 
terms of the density matrix moments.^- A similar expansion is also obtained by computing 
Eq. 0181) explicitly. [The matrix coefficients of the Hamiltonian can be found in appendix 
iBl ] During this computation, it is quite useful to distinguish between atomic kinetic energy, 
Ekin = Tr{(P^/2M)p}, and atomic potential energy, Epot = Trji^eP} (see Eq. ([1])). In terms 
of the CEID moments (see Eq. (!30!l). those two quantities are given by: 



p2 p I 

Ekin = ^ + ^Tr {/io,i} + ^Tre {/io,2} , (33a) 

Epot = Tre |i^e(i?)/io,o} - Tre {F(i?)/ii,o} + ^Tre {kiR)fi2,o} , (33b) 
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As explained in appendix [Cl the CEID evolution of the bare averages defined in Eq. (133|) 
do not give a conserved (bare) total energy, Etot = Ekm+Epot, although the error is negligible 
for large enough CEID order.— That is because CEID provides an approximation of the exact 
evolution (in the truncated Hilbert space) of the observables' averages (see Eqs. IC2l and lC3p . 
On the other hand, the exact evolution (in the truncated Hilbert space) of every observable 
can be retrieved starting from the CEID EOM and then adding a correcting term whose 
general analytical expression is reported at the end of appendix [C] 

The time-derivatives of the corrections for the bare atomic kinetic and potential energy — 
whose integrals must be added to Eq. (l33|) (a) and Eq. (l33l) (b). respectively — are reported 
below: 



CiS„= +^(iV+l)Tr{AF(i?)p^,^} + 

'^°:(iV + 1)^^ Tr [aF{R) {pr,,N-i - Pn-i,n)} + 



AM 



+ 1) V yTr {pr,^N-i + Pn^i,n} + 



{N+l)^TT\^k{R){pN,N-i+PN-i,N)} , (34a) 
+ ^iN+^)^T^^{HR)iPN,N^i-PN^i,N)} . (34b) 



V. COMPARISON BETWEEN EXACT INTEGRATION AND CEID 

In this section we present the main results of this work. They were obtained by means of 
the two numerical algorithms described in the previous sections, namely the exact integration 
scheme of Sec. IIIII and the CEID scheme of Sec. IIVI We recall that our main goal is to 
attest the convergence of CEID (by increasing its order) and to verify that the converged 
results agree with the exact dynamics of the 2LS geometries introduced in Sec. [ITl That 
can be safely done by a direct comparison between CEID and exact integration of the time- 
dependent Schodinger equation and this will be the object of Sec. IV Al — which contains 
a discussion of the electronic observable dynamics — and Sec. IV Bl — which contains a 
discussion of the atomic observable dynamics. 
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The agreement of CEID with exact integration is clearly a fundamental numerical achieve- 
ment, but it does not directly help in the interpretation of the simulation findings. Further 
insights can be obtained by comparing CEID against analytical results derived through first- 
order time-dependent perturbation theory (see appendix [D]). This comparison is reported 
in SecJvn 



A. Electronic observables 



Here we present the results of the electronic dynamics. As initial condition, we always 
choose the atomic vibrational ground-state on the upper PES and than we let the system to 
evolve according to either the exact Schrodinger evolution or the CEID equations. The WT 
of the initial (t = 0) uncorrelated density matrix is: pw{AR, AP,0) = Pq^IAR, AP) pe{0) , 
where Pq^q{AR, AP), according to the definition given in Sec. \IV\ is the WT (in the mobile 
reference frame) of the atomic vibrational ground-state (which is centered in ^ = and 
P = 0) and 

pM = (35) 

describes a pure excited electronic state in the non-adiahatic representation introduced in 
Secini 

The most informative electronic observables are the probabilities to find the system in the 
upper or lower electronic state. Those are obtained as the diagonal entries of the electronic 
density matrix pe = fiQ Q (see Sec. IIV A|) and we shall call them electronic populations. In 
Fig. [T] we collect the numerical results for the unshifted case. 

A sketch of the PES is plotted in the first panel. [In this particular kind of 2LS, the 
difference between adiabatic and non-adiabatic surfaces is negligible.] In Fig. [H^b) the exact 
evolution of electronic populations is reported. We stress that the oscillatory population 
transfer between the two PES clearly confirms the crude picture we guessed in Sec. [Tll 

In Fig. [D^c) the time-evolution of the electronic populations from CEID is reported. 
Different CEID orders, namely = 0, = 1, and A^ = 5, have been considered. As 
we expected, the A^ = (which is equivalent to ED, see Sec. IIV Al) does not display any 
electronic transition. On the other hand, the outcomes of higher order CEID simulations 
present oscillations of the electronic populations which are in almost perfect agreement which 
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FIG. 1: (Color online) Unshifted 2LS. (Top) Adiabatic energy levels of the electronic Hamiltonian, 
Hf,{R) (see Eq. ([2])) are plotted (solid lines) against the atomic coordinate, R. For the units 
employed, see the main text. The horizontal line (red online) marks the the total energy of the 
system. A sketch of the initial atomic density distribution is also given (dashed line). (Center) 
Electronic populations against time, from exact integration. (Bottom) Electronic populations 
against time, from CEID. Results for CEID order (i.e. ED), 1, and 5 are reported. 

the exact integration results. In particular the first order CEID simulation is well converged 
(i.e. there are no visible differences between the = 1 and A^ = 5 findings). On the other 
hand, this is not very surprising; the unshifted 2LS is the easiest case, since the symmetries 
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involved ensure that the evolved state is well described by the simple ansatz stated in Eq. ([3]) 
(see also appendix [Dl) . 

In Fig. [2] we show the results for the shifted case following the same scheme employed for 
the previous figure. For this kind of 2LS, adiabatic and non-adiabatic PES are qualitatively 
different, but since in Fig. [2](a) the difference can be appreciated only close to the crossing, 
we provided a magnified plot of that region in a small inset. Once again, almost perfect 
agreement is seen between a well converged CEID simulation (here for at least = 5) 
and the exact integration of the time-dependent Schrodinger equation, while at the level of 
ED (i.e. = 0) the system is stuck in the upper PES. The fact that a first order CEID 
simulation is not yet well converged is not surprising and is a confirmation of the general 
trend predicted in Sec. [Ill the larger the surface displacement, Rq, the higher will be the 
CEID order required to obtain a well converged simulation. 

We see in Fig. EJ^b) that the period of the electronic oscillations is larger for the shifted case 
than for the unshifted 2LS. [A perturbative account of that effect can be found in appendix 
[Dl ] Moreover, in this shifted case the population exchange between the two PES is not 
complete: the minimum of the electronic population on the upper surface (corresponding to 
the maximum of the electronic population on the lower surface) is not exactly zero (one). 
This interesting feature is clearly visible in Fig. Mjo) and is also found in the two well 
converged CEID simulations in Fig. [21(c) so it is not a numerical feature. This is instead 
a non-trivial fingerprint of otherwise elementary dynamics which is caused by the virtual 
transitions — a clear quantum effect — between the low lying resonant states and more 
energetic atomic vibrational states. Further details can be found in Sec. IV CI in which we 
study the dependence of such residual population on the coupling constant fc- 

B. Atomic observables 

Since the atomic motion is actually non-classical, we expect to find quantum fluctuations 
of the atomic observables around their average values. We study the time-evolution of 
the average atomic position and momentum, R and P, because they provide a sort of 
effective trajectory in the phase space which represent an important link with classical MD. 
Obviously, the concept of trajectory is not well defined in quantum mechanics and is useful 
as an approximation only if the fluctuations are not too large. So, we also consider the 
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FIG. 2: (Color online) Shifted 2LS. (Top) Adiabatic energy levels of the electronic Hamiltonian, 
Hf.{R) (see Eq. ([2|)) are plotted (solid lines) against the atomic coordinate, R. The avoided crossing 
is magnified in the inset, where also the non-adiabatic energies (dashed lines) are reported for 
comparison. For the units employed, see the main text. The horizontal line (red online) marks 
the the total energy of the system. A sketch of the initial atomic density distribution is also given 
(dashed line). (Center) Electronic populations against time, from exact integration. (Bottom) 
Electronic populations against time from CEID. Results for CEID order (i.e. ED), 1, 5, and 10 
are reported. 

variance of the atomic position, (Ai?^), in order to test the accuracy of CEID in describing 
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possibly non-classical atomic dynamics. 



^ 0.05 



-0.05 




-0.1 



FIG. 3: Shifted 2LS. Plots of the time-evolution of the averaged atomic position, R, (top) and 
averaged atomic momentum, P, (bottom). Data (almost perfectly superimposed) are taken from 
exact integration and well-converged CEID simulation. 



We start by reporting results for the average atomic position and momentum, R and P. 
They are evolved by means of the Hamilton-Ehrenfest equations (see Sec. lIVp according to 
CEID while in the exact integration scheme they are obtained by means of Eq. ( ITSl) . For 
the unshifted 2LS, ^ = and P = for all time due to the inversion symmetry displayed 
by the system. On the other hand, the findings reported in Fig. |3] for the shifted case once 
again show almost perfect agreement between CEID and exact integration in a completely 
non-trivial case. In particular, CEID not only reproduce the general trend of both R and P 
(large period oscillations), but also gives the short time scale details (rapid oscillations). 

In Fig. m we report the results for the variance of the atomic position, {AR"^), for both 
the unshifted and shifted 2LS. This observable can been obtained as the trace of the CEID 
moment fi2fl defined in Sec IIV A[ Once again, almost perfect agreement has been found 
between a well converged CEID simulations (here = 5 and = 10 for the unshifted 
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FIG. 4: Plots of the time-evolution of the atomic position variance {AR^), for the unshifted (top) 
and shifted (bottom) 2LS. Data (almost perfectly superimposed) are taken from exact integration 
and well-converged CEID simulation. 

and shifted 2LS, respectively) and the exact integration of the time-dependent Schrodinger 
equation. We stress that such fluctuations are quite significant and so the atomic dynamics 
is only poorly approximated by its average trajectory in the classical phase space. CEID is 
working properly even in those highly non-classical cases. 

Finally, we have also verified the agreement between CEID and exact integration for the 
other entries of the covariance matrix, namely (AP^) and (ARAP). However, numerical 
findings for those cases are nor reported here because they are not qualitatively different 
from the {AR'^) case. 

C. Comparison between time-dependent perturbation theory and CEID 

In this last section we briefly compare the CEID outcomes against time-dependent per- 
turbation theory results. [Mathematical details are collected in appendix [Dl] 
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First of all, from a well converged CEID simulation (e.g. N = 5 and = 10 for the 
unshifted and shifted case, respectively), the values of the electronic oscillation frequency and 
the residual electronic population can be obtained by means of a straightforward numerical 
interpolation. Then, this procedure can be repeated for the same 2LS geometries, but 
different electron-ion coupling constant, fc (see Eq. ([T])). It is instructive to study the 
effect of the atomic motion on the electronic transitions because it might cause non-classical 
phenomena, like quantum interference between different transition paths. Those effects are 
usually hard to interpret without a model which can — at least qualitatively — describe the 
physics involved. Fortunately, for the kind of model 2LS we considered in this paper, a simple 
model can be obtained by means of time-dependent perturbation theory, whose prediction for 
the electronic transition frequency and the residual population are summarized in Eq. flD4p 
and Eq. flD7l) . respectively. 
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FIG. 5: (Color online) Frequency iVp of the electronic transitions (see Eq. (jD4])) (top) and square 
root of the residual population Pris (see Eq. (|D7P ) (bottom) against the coupling constants, fc, for 
the unshifted and shifted 2LS. Linear fits are also showed (dashed hues). For the units employed, 
see the main text. 

In Fig. [5](a) numerical values of electronic population frequencies are reported against 
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several values of the coupling constant, fc- A clear linear trend is manifest in all the 2LS 
geometries. Moreover, numerical values are in almost perfect agreement with the analytical 
results, Eq. flD4D . 

In Fig. \5KJo) the residual populations are plotted against the same coupling constant 
values. Although the analytical trend, P^es — ig"^ (see Eq. (ID7p ) is confirmed, in general is 
not easy to give an estimate of the prefactor, 7. On the other hand, for the unshifted 2LS 
case, only one term in Eq. flD7|l is non-zero due to the SHO selection rules. As a consequence, 
a numerical estimate can be obtained and it gives 7 ~ 2.5 ■ 10~^, while a direct numerical 
interpolation gives 7 ^ 3.7 ■ 10^^. We stress that such disagreement might depend on the 
kind of approximation we made in order to derive Eq. (lD7p and it does not effect the general 
scaling trend of the residual population with the coupling constant. 

VI. DISCUSSION AND CONCLUSIONS 

We have presented a new formulation of correlated electron- ion dynamics (CEID). It 
is based on a suitable expansion of the quantum fluctuations around the mean-field atomic 
trajectories and its lowest accuracy limit has been proved to be equivalent to the well-known 
Ehrenfest dynamics (ED). This new formulation has been obtained by a combined use of: 1) 
an expansion of the density matrix in terms of atomic harmonic states centered around the 
average instantaneous atomic positions; 2) an exact Wigner transform with respect to the 
atomic degrees of freedom of the expanded density matrix. The validity of this scheme has 
been successfully tested by simulating the non-adiabatic time-evolution of a model two level 
system (2LS). The accuracy of our simulations is determined by a single parameter which 
is related to the order of the density matrix expansion and is called the CEID order. We 
then verified that, for all the considered 2LS geometries, the exact quantum dynamics — 
obtained by exact integration of the time-dependent Schrodinger equation — is eventually 
retrieved by increasing the CEID order. We think that this is a crucial property of our new 
CEID scheme which allows us to estimate the convergence of a numerical simulation even 
when reliable benchmarks are not available. 

As for the other proposed CEID schemes,— our algorithm only needs the Hamiltonian 
and the initial conditions to start and the subsequent evolution is computed smoothly, 
without resorting to any kind of surface hopping or wave-function spawning. No a posteriori 
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position, velocity, or density matrix adjustment is needed. The exact evolution (in the 
truncated Hilbert space) of every observable average can be obtained starting from the 
CEID EOM by adding a correction term (whose analytical expression is known) which is 
anyhow negligibly small for large CEID order. Moreover — and at variance with other 
available algorithms^iii — it works perfectly well within a non-adiabatic representation of 
the electronic PES. This is desirable because non-adiabatic PES may be smoother than 
the adiabatic ones^ and also because a costly diagonalization of the atomic potential energy 
Hg{R) at each step is avoided. All these dynamical properties make our new CEID algorithm 
a good candidate for simulating atomic systems in which quantum coherence is relevant. 

The advantages of a coherent quantum scheme might be relevant even when macroscopic 
quantum coherence is not shown. It is well known that ED — at variance with other 
quantum-classical methods^'^- — cannot thermalize a mixed electron-ion system (the elec- 
tronic degrees of freedom are too hot with respect to the atomic ones.'^) This failure of 
the mean-field approximation depends on the absence of quantum fluctuations which cause 
spontaneous phonon emission from an excited electronic state. [This drawback is apparent 
also in the 2LS simulations considered in this paper (see Sec. |Vl): the ED is always stuck in 
the initial excited state.] As a consequence, ED does not satisfy microreversibility.^'^^ On 
the other hand, a CEID simulation beyond ED can describe quantum fluctuations and meet 
the coherence requirements for microreversibility in a very natural way. [Once again, see the 
results of Sec. IVl ] 

Although it is not the main concern of this paper, our group is considering a viable way 
to approach quantum thermalization physics by means of CEID. A first possibility is given 
by the spin-boson model^ in which the bath degrees of freedom are treated explicitly by 
means of a collection of many quantum harmonic oscillators. On a more speculative ground, 
one can think to implement the generalization of the Nose-Hoover thermostat introduced in 



Ref. 



40l . This scheme is known to fail for ED^ due to the lack of correct quantum back- 



reaction on the classical bath variables.— On the other hand, as we have shown again in 
this paper, CEID corrects this ED drawback and it might be better suited for that sort of 
thermostat. Moreover, a successful attempt to couple the Nose-Hoover thermostat to the 
spin-boson model is known in literature^ and it can provide an interesting test case for 
future CEID simulations. 

Our CEID algorithm is computationally demanding and is expected — in the worst case 
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scenario — to scale as (A^+l)^^'=, where is the CEID order and A^^^ is the number of atomic 
coordinates.— Nevertheless, it must be pointed out that the number of relevant atomic 
coordinates can be effectively much smaller than NcM^ In this case, one might accelerate a 
CEID simulation by allowing for quantum atomic fluctuations along the relevant directions 
only. We also stress that the CEID algorithm is still faster than the exact integration scheme 
employed to produce benchmark calculations in this paper (see Sec. IIII|) which should scale 
as (A^ + 1)'^^'= since a numerical diagonalization of the Hamiltonian in the truncated Hilbert 
space is implied. We are also considering alternative truncations of the Hilbert space in 
order to restore the polynomial scaling with atomic degrees of freedom of the early CEID 
algorithms. 

We see another possible advantage of this CEID scheme over exact integration: the former 
expands the quantum fluctuations around mean-field atomic trajectories, while the latter 
expands with respect to a fixed reference frame. Now, consider a quantum motion in which 
there are fluctuations about the mean-field atomic trajectories that are very tightly confined 
along a given direction. With our CEID formulation such fluctuations can be treated ac- 
curately with a low order expansion. However, schemes that employ basis functions which 
are not centered around the atomic trajectories could require a very high order expansion to 
reproduce that confined behavior if the trajectories are remote from the center of the basis 
functions. 

Other algorithms, such as molecular dynamics with quantum transitions or ab initio 
multiple spawning, might have a lower computational complexity, especially if the region of 
the configuration space where non-adiabatic effects are relevant is small and crossed only few 
times during the time-evolution. This is the case, for instance, for many chemical reactions 
in a gaseous or diluted phase. On the other hand, we recall here that CEID was explicitly 
devised to deal with electron-ion correlations in metals, a kind of systems in which the 
aforementioned algorithms are expected to be less efficient. 

Needless to say, a reliable algorithm to simulate microscopic electro-mechanical effects, 
including joule heating, will find important application in nanostructure design. Our CEID 
algorithm is a good candidate because its accuracy can be systematically increased by tun- 
ing a single parameter that allows us to approximate the quantum atomic fluctuations in a 
physically transparent way. Moreover, since quantum coherence is well addressed by CEID, 
subtle photo-physical effect like luminescence in conjugated polymers might be addressed 
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by this method. Apphcations of our algorithm to larger atomic systems and different ther- 
modynamical ensembles are subjects of ongoing study. 
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APPENDIX A: PROPERTIES OF THE P„^(i?,P) FUNCTIONS 

All the properties of the functions Pn,m introduced in Sec. [IV]can be obtained by means 
of the simple harmonic oscillator (SHO) algebra. 
The Pn,m are orthonormal: 

_ 1 f 

{Pn',m'{R, P), Pn,m{R, P)) = T: T / dRdP Pm' ,n' {R, P) Pn,m{R, P) = Sn,n' Sm,m' , (Al) 

zn a J 

and form a complete set in the Wigner space. [Also note that Pn,m = Pm,n-] The following 
identities are useful when computing the action of a canonical operator on Pn,m- 



R Pn,m[Rj P) = +^ 1^/2' '™ y 2 V ^ Rn,m-l + "W ^ Pn,m.+l | ; 

(A2a) 



(A2b) 



pp (pp^- ( hp + i p _ hp I + i p 

n,m l.-f''; -'J 2 I V 2 '™ y 2 2 y 2 | ; 

(A2c) 



/5 p (j^p\- ' I rp -^J'' + ^p hp .h±lp 

'Jp J^n,m\-n", ^ ) — I V "2 '™ y — 2 — "■+1''" ~ y "2" "■''""^ ~ Y — 2 — "■.'^+1 



(A2d) 
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Higher order identities can be recursively obtained by using these basic identities. 



APPENDIX B: DERIVATION OF THE EQUATIONS OF MOTION 

In this appendix we show how Eq. (126|) can be derived from Eq. by using the ex- 
pansion in Eq. ( l25l) . First of all, we want to find the action of the position and momentum 
derivatives on the matrix coefficient Pn,m- We formally introduce the operators D^j and Dp 
in the following way: 

oo oo 

= 5Z 5Z t^"''"] ^"'"^ ' 

71=0 m=0 



oo oo 



^ = E E Dp [Pn,m] Pn,m • (Bib) 
71=0 m=0 

Therefore, by means of the identities reported in Eq. flA2p (b.d). we find that: 



Dfi [Pn,m\ = I Y "2 ~ Y — 2 — ~'~ y "2" ^"''"^^ ~ y — 2 — ^"'"^+1 



(B2a) 



-Up [Pn,m\ ~ [ Y "2 ^""-^'"^ ~'~ y 2 ^^+1'™ ~ \l ~2 ~ \l ^ Pn,m+1 



(B2b) 
(B2c) 

Finally, previous equations can be also written in a more compact form as follow: 

00 

DiJ [pn,m] = J2 [^SPk,ni " Pn,kD^^l , (B3a) 

A:=0 
00 

Dp [pn,rn] = [^2^,,^ - p„,fc^S] , (B3b) 

A:=0 

where Di!% = -{l/ao){>/n/2 6n,m+i - \/m/2 6n+i,m) and Z)i^m = -{i/hQ){^/nj2 5n,m+i + 
5n+i^m)- [Here we employed the Kronecker delta.] As a consequence, the EOM of the 
matrix coefficients pn,m (in the mobile reference frame) can be formally given by a proper 
Lie bracket, as in Eq. (flQll : 

^ oo ^ 

= ^Hn,k Pk,m — Pn,k Hk^m ■ (B4) 



Pn 

'I. a 

k=0 
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where Hn,m = Hn,m + ih{P/M)Di^m +ihFDii^m. [An operator H can be formally defined in 
the original Hilbert space taking the anti-Wigner transform of Yl'^=o X]m=o ^n,m Pn,m-] 

Although Eq. flB4p is completely equivalent to Eq. ([23]), in order to obtain Eq. fl2Bl) we 
still need to compute Hn,m- This can be done by using then properties of the Pn,m{R,P) 
functions (see appendix and the following equation: 

Hn,m = {Pn,m{^R, AP), H^{AR, AP)) . (B5) 

For instance, for a quadratic Hamiltonian expanded with respect to the mobile reference 
frame: 

p2 PAP Ap2 ^ - - - 1 - - 

H^ = + + + He (R) -ARF (R) + -AR^ K (R) , (B6) 

one finds that: 



\/m{m - l)5m-2,n + 



-(2m + l)5m,n + \/n{n - 1)5 

m,n—Z 

^/m{m - l)Sm~2,n + (2m + l)Sm,n + \/ n{n - l)6m,n-2 ■ 

(B7) 

[The first three terms account for the the atomic kinetic energy and the last three for the 
atomic potential energy.] From Eq. ( ]B7I) one can easily obtain Hn,m whose expression must 
be substituted in Eq. (lB4p to obtain Eq. (12^ . The requirement of neglecting in the RHS of 
Eq. fl2Bl) all those matrix coefficients whose indices are greater than the CEID order can be 
directly enforced by constraining the summation index k in Eq. flB4l) to be at most equal to 
the CEID order. 



2M 



M 



m 

2 "m— l,ra 



2 ^m,n—l 



AM 



+ He (R) 6m,n - a^F [R) 



'm~l,n 



+ 



+ ^A- (R) 



APPENDIX C: ORIGIN OF THE CORRECTING TERMS TO THE AVERAGES 



The exact evolution of the average of an observable A is given by the following well-known 
equation: 



A{t) = -Tr <! A 



H,m 



1 

ih 



Tr 



A,H 



m 



(Cl) 
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This expression is still true in the mobile reference frame since the extra terms which arise 
from the implicit time-dependency of A and p cancel out exactly. The truncated version of 
Eq. (ICip is naturally given by: 

1 



ih 



-TriT 



T [pit)] 



(C2) 



where T is the anti-Wigner transform of the truncation operator T^ introduced in Sec. IIV[ 
[We recall that this operator depends on the CEID order.] As expected, averages of those 
observables which commute with the Hamiltonian are constant of motion. 

It turns out that, by using the CEID EOM (see Eq. (!26l) ). the bare dynamics of the 
averages (in the mobile reference frame) is given by: 



in 



+ Tr < T 



T 

dA 

M 



A 



T [pit)] 



P 



T[p(t)]^^ + Tr^T 



OA 



T [pit)] F . 



(C3) 



where we used the modified Hamiltonian H introduced in Sec. [Bl We notice that — unlike 
the situation for the exact evolution — the last two terms in Eq. (ICSP do not cancel out 
exactly the similar terms coming from the modified Hamiltonian H. For this reason, even 
if [T[i], T[/7]] = 0, A'"'"'" is not conserved by the CEID EOM. On the other hand, although 
the Eqs. IC2I and |C3I do not provide the same dynamics, the difference between A^'^\t) and 
^6are^^^ — the Correcting term C^(t) — is expect to be small for large enough CEID order. 
This fact has been verified for the atomic kinetic and potential energy, whose analytical 
expressions of the respective correcting terms are given in Sec. IIVBI 

Finally, we stress that a general expression for the time-derivative of the correcting term 
to the average of the observable A can be derived analytically. It is given by: 



A 



+ Hermitian conjugate , 



(C4) 



where Ot is the Hermitian, idempotent, operator defined by: T[A] = OtAOt- [The operator 
Ot also depends on the CEID order.] By using this expression, one can implement the correct 
truncated dynamics of the observable averages (i.e. Eq. (]C2I) ) starting from the CEID EOM 
and then adding the numerical integral of the appropriate correcting term. 
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APPENDIX D: FIRST ORDER TIME-DEPENDENT PERTURBATION THE- 
ORY OF A RESONANT 2LS 



We start from the non- interacting limit of Eq. ([T]), i.e. by setting = in Eq. ([2]). 
In this easy case, the two PES are dynamically uncoupled. On each surface the atomic 
configurations can be classified by using the appropriate SHO basis set. We indicate by 
IXri*^) and \xn^) the n-th harmonic excitation of the atomic degrees of freedom on the upper 
and lower PES, respectively. In order to obtain the resonance condition between |xo"'') and 
IXi''') described in Sec. [Tll we set their energies equal to Ei = 3/27ja;. [Within this setup, 
the states Ixn""*) and Ixi+i) are also degenerate, having energy En+i = {n + 3/2)huj. The 
ground-state is not degenerate and its energy is Eq = l/2hw] This degeneracy is lifted by 
switching on the interaction (/c > 0). The hybridization energy is: A„+i = {x^n\H\Xj^j^i) ■ 
As a consequence, a system prepared in the state Ixo""*) can decay into the state (see 
Sec. HI])- If the coupling constant is not too large, this process can be treated perturbatively, 
with the small parameter being the adimensional coupling constant g = fcdo/huj. 

We define the unperturbed Hamiltonian Hq by neglecting all the entries of H except the 
diagonal energies and the two off-diagonal matrix elements coupling Ixo""*) and 



Ha 



Eo ■■■ 

El Ai 

Ai El 

: ^2 



V 



:dii 



The perturbation matrix given hj V = H — Hq has zero diagonal entries and does not couple 
IXo"^) and \x?). 

According to our initial condition, the zero order solution of the time- dependent 
Schrodinger equation ihdt\'^{t)) = ^|\E'(t)) is given by: 



cos 



|Ai|t 



Xo ') - «sm 



|Ai|t 



I {0\ 
IXi ) 



From Eq. flD2p . it is easy to derive the dynamics of the upper electronic population: 



cos 



|Ai|t 



(D2) 



(D3) 
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The correspondent oscillation frequency is given by: 

Up = 2-^-^ (X gu (D4) 

We stress that at half period {t = T/2 = n/up) the upper electronic population is exactly 
zero, i.e. at zero order there is no residual population (see Sec. IV Ap . 

To go beyond the zero order, we use the following transformation: \'^{t)) = e^^°*\(j){t)) 
which leads to an effective Schrodinger equation for 10): 

ihdtm)) = vitmt)) , (D5) 

where V{t) = e"«ft^«*V^e«^o*. This is an interaction picture transformation. Eq. (IDSP can 
be solved using the standard iterative procedure (Volterra's equation) giving the following 
first order approximation of the wave-function: 

\¥'\t)) = \¥'\t)) + ^ dr V{t)\xP) . (D6) 

The last term on the RHS of Eq. flD6p generates the finite residual population at half 
period observed in our numerical experiments. The following simplified expression holds if 
hu > |Ai|: 

I / I T/ 1 O \ 1 2 

P.e. = P«(T/2) = I (xf^|vI/«(T/2))p = J] ' '^'^i oc . (D7) 
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